arXiv: 1503.05319vl [physics.flu-dyn] 18 Mar 2015 


Non-linear Langevin model for the early-stage dynamics of 

electrospinning jets. 

Marco Lauricella^, Giuseppe Pontrelli^, Dario Pisignano^’^, and Sauro Succi 

^Istituto per le Applicazioni del Calcolo CNR, Via dei Taurini 19, 00185 Rome, Italy 
^Dipartimento di Matematica e Fisica Ennio De Giorgi,University of Salento, via Arnesano, 

73100 Lecce, Italy 

^Istituto Nanoscienze CNR, Via Arnesano 16, 73100 Lecce, Italy 
Wednesday 28*'^ January, 2015 


Abstract 

We present a non-linear Langevin model to investigate the early-stage dynamics of electrified polymer 
jets in electrospinning experiments. In particular, we study the effects of air drag force on the uniaxial 
elongation of the charged jet, right after ejection from the nozzle. Numerical simulations show that the 
elongation of the jet filament close to the injection point is significantly affected by the non-linear drag 
exerted by the surrounding air. These result provide useful insights for the optimal design of current and 
future electrospinning experiments. 


1 Introduction 

The dynamics of charged fluids under the effect of an external electrostatic fields represents a major theme of 
non-equilibrium thermodynamics and statistical mechanics [I1121E1I1]. 

Charged liquid jets may develop several types of instabilities, depending on the relative strength of the various 
forces acting upon them, primarily electrostatic Coulomb self-repulsion, viscoelastic drag, surface tension effects 
and dissipative forces due to the interaction with the surrounding environment. Since these instabilities are 
central to many industrial processes, including the production of ultrathin fibers via the so-called electrospinning, 
they have been analyzed since the late 1960s [S]. However, due to inherent difficulties associated with the 
underlying long-range many-body problem, a thorough theoretical understanding is still lacking, whence an 
important role for computer simulation. 

In the recent years, the production of ultrathin nanofibers has found increasing applications in micro and 
nanoengineering and life sciences as well [a m. In the electrospinning process, electrospun nanofibers are 
produced at laboratory scale by the uniaxial stretching of a jet, which is ejected at the nozzle from the surface of a 
charged polymer solution. This initial elongation of a jet is produced by applying an externally electrostatic field. 
An intense electric field (typically 10^ — 10®V -m"^) is applied between the spinneret and an oppositely charged 
collector, usually placed at about 20 cm from the injector. Electrospinning involves mainly two sequential 
stages in the uniaxial elongation of the extruded polymer jet: an initial growth stage, in which the electric 
field stretches the jet along a straight path away from the nozzle of the ejecting apparatus, and a second stage 
characterized by a bending instability induced by small perturbations, which misalign the jet away from its 
axis of elongation. These small disturbances may originate from Coulomb repulsion on different portions of the 
jet, as well as from mechanical vibrations at the nozzle or aerodynamic perturbations within the experimental 
apparatus. As a consequence, the jet path length between the nozzle and the collector increases, and the stream 
cross-section undergoes a further decrease. The prime goal of electrospinning experiments is to minimize the 
radius of the collected fibers. By a simple argument of mass conservation, this is tantamount to maximizing the 
jet length by the time it reaches the collecting plane. Consequently, the bending instability is a desirable effect, 
as it produces a higher surface-area-to-volume ratio of the jet, which is transferred to the resulting nanofibers. 
By the same argument, it is therefore of interest to minimize the length of the initial stable jet region. 

Simulation models provide a useful tool to elucidate the phenomenon and provide valuable information for 
the design of future electrospinning experiments. Numerical simulations enhance the capability of predicting 
the key process parameters and exert a better control on the resulting nanofiber structure. In recent years, with 
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a renewed interest in nanotechnology, electrospinning studies attracted the attention of many researchers both 
from modeling and experimental points of view UM- These models treat the jet hlament either as a charged 
continuum fluid, or as a series of discrete elements (beads) obeying the equations of Newtonian mechanics[7]. The 
latter standpoint is the one taken in this work. Each bead is subject to different types of interactions, namely 
long-range Coulomb repulsion, viscoelastic drag and the external electric field. The main aim of such models 
is to investigate the complexity of the resulting dynamics and provide the optimal set of parameters driving 
the process. Actually, due to the large number of experimental parameters, electrihed jets are still treated via 
empirical approaches. The effect of fast-oscillating loads on the bending instability, have been explored in a 
modeling and computational study m- On the other hand, the effect of air drag in electrospinning process has 
been addressed only very recently HU, even though there is experimental evidence that the air drag affects the 
dynamics of the nanohber via a non-linear dependence on the jet geometry [12]. 

Here, we investigate the uniaxial elongation of an electrihed polymer jet in the early-stage of its dynamics at 
the nozzle of the ejecting apparatus and under stochastic-dissipative force. This perturbation effect is modeled 
by a Langevin approach. In particular, we assume that a Brownian term can adequately reproduce a stationary 
perturbation related to multiple simultaneous tiny impacts along the direction of the jet elongation axis, as in the 
case of an air drag force generated by the motion of a polymer jet through a gaseous medium. Relations between 
air drag force and Brownian motion were proposed in literature starting from experimental observations [13lll4j . 
In a recent work, we extended the one-dimensional bead-spring model developed by Pontrelli et a/.|15j. to 
include a linear dissipative-perturbing Langevin force which models the effects of the air drag force. This is 
accomplished by adding a random and a dissipative force to the equations of motion, and further assuming 
that the two terms obey the huctuation-dissipation theorem m- As a result, the system was described by a 
Langevin-like linear stochastic differential equation HU- 

Based on empirical evidences, in this paper we extend the previous approach to include non-linear effects 
due to air disturbances (see Sec 2). The resulting non-linear Langevin equation is numerically integrated to 
investigate different dynamical regimes under the effect of systematic and random air perturbations (see Sec 3). 


2 The mathematical model 


Let us consider a rectilinear electrihed viscoelastic jet in a typical electrospinning experiment. In order to model 
the stretching, we represent the hlament by a viscoelastic dumbbell (or dimer) with two charged beads (a, b) 
with the same mass m and charge e (not to be confused with the electron charge) located a distance I apart 
(Fig. g. One of the two beads (a) is held fixed to the nozzle, while the other one (b) is free to move under 
the effect of the various forces. The dumbbell ab represents a schematic model for the jet. We denote by h 
the distance of the collector plate from the injection point and by Vq the applied voltage between them. The 
bead b is subject to three different forces (gravity and surface tension are neglected): the external electrical 
field Vo/h, the Coulomb repulsive force between the two beads, and the viscoelastic force. As pointed out 
elsewhere HU 13 El, the competition between Coulomb and viscoelastic forces characterizes the early-stage of 
the jet elongation, while the second stage is dominated by the external electrical force, driving the fluid hlament 
towards the collector. 

The combined action of these three forces governs the elongation of the jet according to the following 
equation jH]: 


dv eVo n 


( 1 ) 


where v is the velocity of the bead b, t the time, r the cross-sectional radius of the hlament, Trr^cr the force 
pulling the bead b back to a, with a the stress of the viscoelastic force. Assuming a viscoelastic Maxwellian 
liquid [7] the time evolution of the stress a related to the viscoelastic force is provided by the equation: 


da ^dl G 
at Idt fa 


( 2 ) 


where G is the elastic modulus, and // the viscosity of the huid. The velocity v satishes the kinematic relation: 


(U 

dt 


= V. 


(3) 


To include air drag effects on the dynamics of electrihed jets, we add a random and dissipative term into 
Eq (H|. Denoted Dy the diffusion coefficient in velocity space and a the friction parameter, we assume that the 
dissipative term has the form aW, while the random force term the form y/OU^rj (t), where ry (t) is a nowhere 
differentiable stochastic process with < ry (H) ry (^ 2 ) >= d {\t 2 — D|), and < ry (t) >= 0. As a hrst approximation, 
we have taken a constant HU- Adding these two force terms in Eq Q, we obtain a Langevin-like stochastic 
differential equation: 
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dv eVo n 

m— = TTT H—;- Trr a — mav 

dt r h 


+ \/2mPlD^rj (t) 


(4) 


A dependence of the dissipative term both on the velocity and on the length of nanohbers was experimentally 
reported czmiiz], and the following empirical formula for the dissipative air drag force was proposed: 


fair = (5) 

where a denotes an empirical factor depending on the air density and kinematic viscosity of gaseous medium. 
Therefore, we propose here a more comprehensive model for the Langevin-like stochastic differential equation 
([^3) which takes the form: 

= — 7rr^(7 — mal {t) , (6) 

dt n 

where the parameter P denotes the super-linearity of the dissipation term, and S accounts for a sub-linear 
dependence of the dissipative term on the jet length. It is worth observing that Eq ([^ reduces to Eq Q in the 
limit P,S^0. 

All the variables and equations are recast in a more convenient non-dimensional form [H]. To this aim, we 
define a length scale L = [e ^, at which Coulomb repulsion matches the reference viscoelastic stress 
G, being ro the initial radius, and r = /r/G a relaxation time. By defining 
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and applying the conservation of the jet volume, vrr^ I = the equations of motion in terms of non- 

dimensional (barred) variables take the following form: 
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2D^7] (i) 


(8a) 

(8b) 

(8c) 


The dimensionless groups are given by: 


"^0^' (9) 

hLmG^ 

2 2 

mLG 

The three group measure the strength of Coulomb interactions, external field and viscoleastic forces in units 
of the inertial force mdv/dt, respectively. In these units, F^e. = Qj so that we are left with four independent 
parameters (considering also a and D^). 


Q = 

V = 


2.1 Numerical integration 

The equations of motion equations are discretised on a uniform sequence ti = tg + iAt, i = 1,... ,nsteps- At 
each time step, we first integrate the stochastic Eq ([^) using the explicit “strong” order scheme proposed by 
E. Platen [T71118] . whereof the order of “strong” convergence was evaluated in literature equal to 1.5. 

The Platen scheme delivers: 
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( 10 ) 


with denoting the approximation for the fc-th component of a generic vector X whereof the time derivative 
is dX/dt = a (t, ..., + b dVt, denoting ft (t) a Wiener process. Here, the vector supporting values T± 

are 


f± = Yt + aAt±b(At)'/% (11) 

and Aft and A*P' are normally distributed random variables related to two independent N (0,1) standard 
Gaussian distributed random variables Ui and U 2 via the linear transformation: 


Aft = Ui (At)'/" 



( 12 ) 


Note that for 


b 


0 the Platen scheme reduces to the second-order Runge-Kutta scheme. Finally, the 


remaining Eqs ([^) and ([^) are integrated via a second order Runge-Kutta scheme with the Hi+i value 
previously obtained via the Platen scheme. 


3 Results and Discussion 

We investigate different regimes of electrihed jets associated with the Eqs ([^, with special focus on metastable 
states and asymptotic behavior. As a reference case, we consider the typical values of Q = 12 and V = 2, 
already investigated in previous works [TKl |S]. All simulations start from the same initial conditions: 1 = 1, 
d = 0 and IK = 0. Note that for reference parameters developed by experimental resultsjH] the typical values 
of length scale L and relaxation time r are 3.19 mm, and 10“^ s, respectively. 


3.1 Parameter setup and asymptotes 


First, we study the deterministic case, by imposing d = = 0. We integrate in time forward and backward 

Eqs in the interval G = 0 and 4 = 5 at different values of time-step At in order to assess a suitable value for 
the specihc case under investigation. Exploiting the time reversibility, we measured an average absolute error 
Al = Ihristeps ~ ^o| lower than 10“^^ with time step At = 10“^. Therefore, we take a time step At = 10“^, as 
a conservative choice for the specific case under investigation. 

We next discuss the elongation of the jet under stochastic perturbation. To estimate the parameters d, S 
and P, we consider the empirical formula Eq ([^ for the dissipative air drag force [121 [3 E]. For typical air 
density, kinematic viscosity of air gaseous medium and jet mass (assumed constant), we obtain a d ~ 0.5 (with 
T = 10“^ s), while the parameters S and P are set equal to 0.905 and 0.19, respectively. For all the investigated 
cases, we adopt for the sake of simplicity the value = d. In order to accumulate sufficient statistics, we have 
run 10000 independent trajectories for each different case under investigation. Thence, we compute the time 
dependent mean value of our observables along the dynamics. 


We integrate the Eqs (8c I for three different cases: in the first, we set d = 0, 5 = 0 and P = 0 {deterministic 
case); in the second d = 0.5, 5 = 0 and P = 0 {linear Langevin); for the third case {non-linear Langevin) 
d = 0.5 while and 5 = 0.905 and P = 0.19, respectively. 

Three basic time-asymptotic regimes can be identihed: 

i) Deterministic: W (x P, I (x P. This is the free-fall regime driven by the external voltage, once every other 
force is extinguished (see Fig[^. 

ii) Stochastic, Linear Dissipation: W oc const, I oc t. This is the ballistic regime resulting from the balance 
between the external held and linear dissipation. Asymptotically, the jet moves at a constant speed, like electrons 
in a linear host media. 

iii) Stochastic, Non-linear: I ex PP, W oc t~^P. This is the regime resulting from the balance between the 
external held and non-linear dissipation with 5 = 0.905 and P = 0.19. Note that, even though these exponents 
are close to the case of standard diffusion, they result from a very different process, namely a constant force 
against a drag growing with the hlament length. It is interesting to notice that the hlament length is still 
unbounded in the limit t —>■ 00 (see Fig |^, but much slower than the ballistic case associated with linear 
dissipation. We emphasize that these asymptotic regimes, although important for a qualitative analysis of the 
process, bear limited practical interest since, in the long-term, the jet undergoes a bending instability which 
cannot be described by the present one-dimensional model. The one-dimensional model is however very useful 
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to discuss the early-stage of the evolution and incipient onset of the instability. In Fig|^ we report the velocity 
W (f) versus time for the three cases above. 


3.2 Deterministic case: no dissipation 

In the deterministic case, we identify two sequential stages in the elongation process (denoted A and B in 
Fig 3^. In the hrst regime, we observe a small increase of W (t^, which rises up to achieve a quasi stationary 


point denoted by t*, where the viscoelastic force 




balances the sum of the two force terms 


Q 


^ i^*) I (t*) 

and V, providing a zero total force. Subsequently, after about 20 — 40 ms, the velocity attains a near-linearly 
increasing trend, close to the time-asymptotic solution discussed above (see also Ref m)- Note that the instant 
corresponds to the lower limit of the derivative dW (Tj /dt, and discerns the two stages of the dynamics, 
regime A characterized by the competition between Coulomb and viscoelasticity, and regime B, characterized 
by the sole action of the external held. 


3.3 Stochastic case: linear dissipation 

Once linear dissipation is included, the jet elongation suffers an additional slow down, leading to a decrease 
of the velocity after the early peak driven by Coulomb forces. This leads to a local minimum at t** with no 
counterpart in the deterministic case, as already discussed in Ref m- We observe that the time occurrence 
of the peak is anticipated by the air drag, an the corresponding velocity slightly decreased. Subsequently, in 
the regime B the velocity attains its asymptotic value (see Fig§. 

3.4 Stochastic case: non-linear dissipation 

In the third case (non-linear Langevin), we observe that the non-linearity of the dissipative term largely alters 
the time evolution of velocity. In particular, it leads to just one quasi stationary point. Furthermore, the 
velocity W (t) appears to tend asymptotically to zero, as it pertains to its asymptotic regime 

Q 


We next examine the force terms, as shown in Fig 


First, the Coulombic term 




decays rapidly so 


that it plays no role in regime B. The early-stage of the dynamics is characterized by the terms 


FveO' {i) 


and 


al , which increase as consequence of the jet stretching due to the external electric held. This early-stage 


comes to the quasi stationary point at time t*, where 


FveO- (t) 

uj) 


and conspire to balance the external 


drive V. 

After the instant f*, the term becomes lamer in magnitude, leading to a decrease of the velocity 


IT, due to the dissipative and viscoelastic terms (see Fig 


3). At the same time, the term —=■ 


Fvefy (i) 


Jit) 


starts to decay. 


and the jet dynamics is governed only by the remaining opposite terms V and . Since dissipation 

grows quasi-linearly with the jet elongation, the velocity goes to zero in the time asymptotic limit. 

Next, we investigate the elongation of jet under stochastic perturbation modelled by the non-linear Langevin 
equation for different values of a, keeping the condition = a. In particular, we explore the way that the 
position of is altered by the dissipative-perturbing force + \/21)^1] (F) in Eq (8cl. In Fig[^ 


we 


report the time evolution of the velocity W (F) for different a. For all investigated regimes we observe a quasi 
stationary point at time F*, which decreases by increasing the term a (see TabiR. 

The straight path of the electrihed jet is described by the observable F (F*j7 which is seen to decrease by 
increasing a. (Fig[^. Useful hints for the optimal design of the electrospinning processes, resulting from deeper 
insights into the early-stage dynamics of the jet can be numerous, including the possibility of better controlling 
the subsequent development of three-dimensional instabilities, and consequently, the diameter and morphology 
of collected nanostructures, as well as the assembly and positions of nanohbers impinging onto the collector. 
A better control of these nanofabrication processes could therefore entail the identihcation and tailoring of air 
drag mechanisms, which can eventually be induced by a proper designed of gas-injecting systems nearby the 
spinneret. 


4 Conclusions 

Summarizing, we have investigated the flow of charged viscoelastic fluids in the presence of stationary stochastic 
perturbations. A Brownian term has been used to model the effects of a perturbing force on the stretching 
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properties of electrically charged jets, providing significant qualitative new insights. We demonstrated that the 
non-linear dependence of the dissipative term on the geometry of the electrospun polymer provides important 
effects that cannot be properly modelled by a linear Langevin-like stochastic differential equation. We also 
observed that the air drag force significantly affects the dynamics of the electrospinning process, leading to 
a time-asymptotic vanishing velocity of the jet. Furthermore, a reduction of the linear extension of the jet 
is observed in the early-stage by increasing the dissipative force term. These results may contribute to the 
optimal set-up of the experimental conditions, so as to enhance the efficiency of the process and the quality of 
the electrospun hbers. These may include, among others, environmental vibrations and resulting micro-vorticity 
patterns. 
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Tables 


a 


1 (t*) 

d (t*) 

w{u) 

0 

0.86 

3.39 

0.81 

3.52 

0.01 

0.81 

3.21 

0.80 

3.46 

0.05 

0.64 

2.58 

0.71 

3.32 

0.1 

0.55 

2.29 

0.65 

3.20 

0.5 

0.40 

1.75 

0.47 

2.69 


Table 1: Values of the dimensionless variables length I, stress a and velocity W at the quasi-stationary points 
computed for different values of a. To be noted the decrease of at increasing a .Furthermore, we stress 
the decrease of the lenght I (t^) (see Fig. [^, proving that the initial stage of the elongation process contracts 
as a consequence of the uniaxial perturbation. 


Figures 



Figure 1: Schematic drawing of the electrospinning process (not in scale), showing h the distance between the collector 
plate and the injection point (nozzle), Vo the applied voltage between these two elements, and the z reference axis whose 
origin is fixed at the injection point. 
























Figure 2: Time evolution of the jet elongation I (t) for three different cases: 1) Deterministic: d = 0, S' = 0 and P = 0 
in black continuous line; 2) Linear Langevin: a = 0.5, S = 0 and P = 0 in red dashed line; 3) Non-linear Langevin: 
a = 0.5, S = 0.905 and P = 0.19 in blue dotted-dashed line. In green dotted line we report the linear fitting (in log-log 
scale) for the Non-linear Langevin case with slope equal to the expected value ‘^/r. The horizontal line on the left side 
corresponds to the physical length 1 = 3.19 cm. 



Figure 3: Time evolution of the velocity W (tj for three different cases: 1) Deterministic: d = 0, S = 0 and P = 0 
in black continuous line; 2) Linear Langevin: d = 0.5, S = 0 and P = 0 in red dashed line; 3) Non-linear Langevin: 
d = 0.5, S = 0.905 and P = 0.19 in blue dotted-dashed line. Two stages of the elongation process are identified: an early 
transient (A), which comes to a quasi stationary point (denoted by a star symbol) and a later stage (P) controlled by 
the competition between the external field and air dissipation. The figure highlights the major role played by dissipation 
on the long-term evolution of the system, which is from free-fall acceleration to a constant velocity regime. As expected, 
the super-linear dissipation enhances the drag effect, leading to a further reduction of the time-asymptotic speed. The 
horizontal line on the right side corresponds to the physical speed v = 79.9 cm/s. 
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Figure 4: Time evolution of the four force terms versus time t: (continuous line), (dashed line), V 


Q 


i 


(dotted line), and — (dashed-dotted line) for the non-linear Langevin case a = 0.5. The quasi-stationary point t* = 0.4 

_ 

is highlighted. The viscoelastic force peaks at about t = 0.5. Subsequently, the force terms a,nd ^ decay 

to zero, while the term and the external electrical field V govern the jet dynamics. 



Figure 5: Time evolution of the velocity W (t) for different values of d. From top to bottom curves: d = 0 (continuous 
line), 0.01 (dashed line), 0.05 (dashed-dotted line), and 0.1 (dotted line), and 0.5 (dashed-dotted-dotted line), keeping 
Di, = d for all the cases. The quasi-stationary points are denoted by a star symbol. As expected, increasing the drag 
coefficient entails a substantial reduction of the filament velocity. The horizontal line on the right side corresponds to 
the physical speed v = 127.6 cm/s. 



Figure 6: Time evolution of the jet elongation I (tj for different values of d: 0 (continuous line), 0.01 (dashed line), 0.05 
(dashed-dotted line), and 0.1 (dotted line), and 0.5 (dashed-dotted-dotted line), keeping = d for all the cases. The 
quasi-stationary points are depicted as star symbol. The horizontal line on the right side corresponds to the physical 
length I = 0.638 cm. 
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Supplemental Information 

In this paper, we use a recast form of the empirical formula for the air drag force fair usually found in literatur^ 
Here, we provide few details on our revised expression. The original empirical formula for the air drag force 

fair is 


(2vr\ 

fair = I ■ O.GbTrrpairv'^ ( - ] , ( 1 ) 

where we consider I as the distance between the beads labeled i and z — 1 , r the cross-sectional radius of the 
filament, v is the velocity of the i — th bead, pa is the air density, and I'a is the air kinematic viscosity. 

We rearrange the Eq as 


/ 2 \ 

fair = I ■ O.eSvrpair I - 1 (2) 

\ ^air J 

Assuming a constant volume of the jet = TTr^L, so that r = with L and tq respectively the 

length and the radius of the jet segment between the beads i and z — 1 at the nozzle before the stretching, we 
obtain 


Thus, we can write 


If we write fair as 


/ 2 

fair = I -j I 

\ ^air J 


0.57-0.5\0-19 _i,19 


fair 


O.Gbnpair 1 - 

^n.i.r 


-0.81 


rO.095.^0.19 
Lj Tn 


0.905„,1.19 




(3) 

(4) 




(5) 


denoting rrii the mass of the bead z, we obtain an equivalent formula for the air drag force with the dissipative 
coefficient a equal to 


a = O.eSvrpair 


2 

^air 


-0.81 


LO 


.095m.19 


( 6 ) 


^A. Ziabicki and H. Kawai, "High-Speed Fiber Spinning: Science and Engineering Aspects", Krieger Publishing Co (1991); 
Sinha-Ray, Suman and Yarin, Alexander L and Pourdeyhimi, Behnam, "Meltblowing: I-basic physical mechanisms and threadline 
model", Journal of Applied Physics (2010), 034912; Yarin, Alexander L and Pourdeyhimi, Behnam and Ramakrishna, Seeram, 
"Fundamentals and Applications of Micro and Nanofibers", Cambridge University Press (2014). 
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